Fast estimation of adult cerebral blood content and oxygenation with hyperspectral time-resolved near-infrared spectroscopy

Near-infrared spectroscopy (NIRS) can measure tissue blood content and oxygenation; however, its use for adult neuromonitoring is challenging due to significant contamination from their thick extracerebral layers (ECL; primarily scalp and skull). This report presents a fast method for accurate estimation of adult cerebral blood content and oxygenation from hyperspectral time resolved NIRS (trNIRS) data. A two-phase fitting method, based on a two-layer head model (ECL and brain), was developed. Phase 1 uses spectral constraints to accurately estimate the baseline blood content and oxygenation in both layers, which are then used by Phase 2 to correct for the ECL contamination of the late-arriving photons. The method was validated with in silico data from Monte-Carlo simulations of hyperspectral trNIRS in a realistic model of the adult head obtained from a high-resolution MRI. Phase 1 recovered cerebral blood oxygenation and total hemoglobin with an accuracy of 2.7 ± 2.5 and 2.8 ± 1.8%, respectively, with unknown ECL thickness, and 1.5 ± 1.4 and 1.7 ± 1.1% when the ECL thickness was known. Phase 2 recovered these parameters with an accuracy of 1.5 ± 1.5 and 3.1 ± 0.9%, respectively. Future work will include further validation in tissue-mimicking phantoms with various top layer thicknesses and in a pig model of the adult head before human applications.


Introduction
Near-infrared spectroscopy (NIRS) is a portable technology that uses safe (i.e., non-ionizing) near-infrared light to noninvasively probe living tissue (Denault et al., 2018). NIRS has high sensitivity to key biomarkers of brain health such as cerebral blood content and oxygenation (Durduran et al., 2010), and is now widely used for neuromonitoring in both pre-clinical and clinical settings (Brown et al., 2002;Goldman et al., 2004;Lewis et al., 2018;Abdalmalak et al., 2020a,b;Milej et al., 2020b). Several NIRS methods have been developed over the years, but the most popular are based on continuous-wave NIRS (cwNIRS), which is the simplest NIRS technology and is based on sending a light beam of constant intensity into the tissue Frontiers in Neuroscience 01 frontiersin.org and monitoring for changes in light attenuation (Scholkmann et al., 2014;Wojtkiewicz et al., 2014). Importantly, by using a cwNIRS technique that can measure light attenuation at dozens of wavelengths, rather than just a few, the spectral features of the main light absorbers in tissue (e.g., water, oxyhemoglobin and deoxyhemoglobin) can be leveraged to improve the accuracy of cwNIRS (Diop et al., 2014). This approach is often called hyperspectral cwNIRS and has been shown to be reliable for neonatal neuromonitoring; (Diop et al., 2014;Rajaram et al., 2018) however, hyperspectral cwNIRS neuromonitoring in adults remains a challenge due to significant contamination from their thicker extracerebral layers (ECL; scalp, skull, and cerebrospinal fluid) Milej et al., 2020a;Li and Diop, 2022).
To mitigate this challenge, alternative NIRS methods that are more sensitive to deep-lying tissue have been developed, and the most advanced are based on time-resolved NIRS (trNIRS) (Gagnon et al., 2008;Diop et al., 2010). In trNIRS neuromonitoring, short pulses of light are released into the head, and the arrival time of each photon at the detector is precisely measured to generate a distribution of timeof-flight (DTOF). This allows for the differentiation of early arriving photons, which have only passed through the ECL, from late-arriving photons which are more likely to probe the brain. Several studies have demonstrated the superior brain sensitivity of trNIRS; (Chance et al., 1988;Liebert et al., 2003;Selb et al., 2005;Steinbrink et al., 2006) however, analyzing trNIRS brain measurements is challenging. Consequently, measurements are often interpreted by assuming that the brain is a homogeneous medium (Patterson et al., 1989;Gerega et al., 2018;Baker et al., 2019). While such an approach may work in neonates, given their thin (3-5 mm) ECL (Li et al., 2015;Sharma et al., 2020), modeling the adult head as a homogeneous medium is an oversimplification since NIRS probes are typically positioned on the scalp for non-invasive measurements. As such, adult measurements contain significant contribution from the ECL since light must travel through skin (∼1.2-1.5 mm) and thick skull (∼10 mm) before reaching the brain (Selb et al., 2014b).
A simple, yet more accurate, approach to account for the ECL contribution is to divide the head into two compartments: brain and ECL (Gagnon et al., 2008). This permits the use of analytical solutions of light propagation in two-layer turbid media to analyze adult trNIRS brain measurements (Kienle et al., 1998). Nevertheless, analyzing trNIRS data with a two-layer analytical model is prone to crosstalk because of the increased number of fitting parameters (from 3 for a homogeneous medium, to 6): 4 parameters for the absorption and scattering coefficients of both layers, the thickness of the top layer, and an amplitude term that accounts for the unknown gain of the trNIRS system. We hypothesize that trNIRS data at dozens of wavelengths, similar to hyperspectral cwNIRS, will allow for the use of the spectral features of tissue chromophores to better constrain the fitting and reduce crosstalk.
To test this hypothesis, we developed a two-phase fitting algorithm−based on a two-layer analytical model of light propagation in diffuse media (Kienle et al., 1998)−that leverages the spectral features of oxyhemoglobin and deoxyhemoglobin to estimate their concentration in the adult brain from hyperspectral trNIRS data acquired with probes positioned on the scalp. Phase 1 of the algorithm uses the two-layer analytical model to accurately estimate initial (i.e., baseline) chromophore concentrations, and Phase 2 uses these initial values and late-photon analysis (i.e., the tail of the DTOFs) to rapidly estimate subsequent concentrations. This approach significantly reduces computation time while accounting for the ECL contribution to the optical signal. The accuracy of the method was validated with in silico data from Monte-Carlo simulations of hyperspectral trNIRS in a realistic model of the adult head obtained from a high-resolution MRI. The estimated concentrations were compared with the known inputted values using Pearson's correlation and Bland-Altman plots.

Materials and methods
The different steps of the two-phase fitting algorithm are illustrated in Figure 1. In Phase 1, data from two source-detector distances (2 and 3 cm) are fit separately with a solution to the diffusion approximation (DA) for a semi-infinite two-layer medium to recover the absolute concentrations of oxyhemoglobin and deoxyhemoglobin in the brain (i.e., baseline concentrations) (Kienle et al., 1998). Phase 2 focuses on rapid estimation of the brain chromophore concentrations once their baseline values are known. It is noteworthy that Phase 2 uses data from both source-detector distances to account for the ECL contribution to the signal.

Phase 1: Estimation of baseline chromophore concentrations
There are typically six fitting parameters when a trNIRS curve (i.e., DTOF) is analyzed with an analytical model of light transport in a two-layer diffuse medium: the scattering and absorption coefficients of each layer, the top layer thickness, and an amplitude term. Because of the large number of fitting parameters, such a procedure is prone to crosstalk between the parameters (Gagnon et al., 2008). The current algorithm reduces crosstalk by dividing the fitting procedure into five steps: (i) roughly estimate the ECL chromophore concentrations, (ii) roughly estimate the brain chromophore concentrations, (iii) estimate the scattering coefficients, (iv) refine the ECL concentrations, and (v) refine the brain concentrations. Further, the algorithm leverages the spectral content of the data to better constrain the fitting, thereby further reducing crosstalk. The details of each step are provided in the following sections. 2.1.1.
Step 1: Rough estimation of the concentrations of oxyhemoglobin and deoxyhemoglobin in the ECL The goal of this step is to use the shorter (2 cm) source-detector distance data to obtain rough estimates of the concentrations of oxyhemoglobin and deoxyhemoglobin in the ECL, which will be further refined in Step 4. Since more than 80% of the signal obtained at the 2 cm source-detector distance comes from the ECL , analyzing this data with a solution to the DA for a semi-infinite homogeneous medium will yield results that are heavily weighted toward the ECL parameters. The homogeneous medium fitting was implemented in MATLAB 2020B using a bounded least-squares regression method [fminsearchbnd (D'Errico, 2021)] to minimize the difference between the analytical model of light propagation in a semi-infinite homogeneous medium and the DTOF data at each wavelength. The DTOFs were fit from 50% of the max on the leading edge to 5% of the max after the peak. The three fitting parameters in the homogeneous fitting were the absorption Flowchart of the two-phase hyperspectral trNIRS data analysis algorithm. The following provides additional information on each step of the algorithm.
Step 2: The tail of the logarithm of the DTOFs was fit to obtain µ a (λ); Computed [Hb] and [HbO 2 ] from µ a (λ).
Step 4: [Hb] and [HbO 2 ] in the brain were fixed to the values obtained in Step 2; µ s '(λ) fixed to values obtained in Step 3; Fitting parameters: ECL thickness, and [Hb] and [HbO 2 ] in the ECL; The initial guess of [Hb] and [HbO 2 ] were the values obtained in Step 1 and were allowed to vary ± 50%; ECL thickness was set to the average value obtained from the segmented MRI and allowed to vary ± 10%.
Step 5: [Hb] and [HbO 2 ] in the ECL were fixed to the values obtained in Step 4; µ s '(λ) fixed to values obtained in Step 3; Fitting parameters: [Hb] and [HbO 2 ] in the brain; The initial guess of [Hb] and [HbO 2 ] in the brain was 30 µM and allowed to vary between 0 and 80 µM.
Step 6: The absorption coefficients in the ECL, µ a ECL (λ), and the brain, µ a Brain (λ), obtained from Steps 4 and 5 were used to estimate the ECL contribution (f ECL ) to the absorption coefficients from the late-photon analysis, µ a Lp (λ).
Step 7: For any subsequent trNIRS spectrum (i.e., DTOFs dataset), the absorption coefficient in the brain, µ a Brain (T, λ), can be rapidly estimated from the baseline absorption coefficients in the brain and ECL, and the ECL contribution to the absorption coefficients obtained from late photon analysis of the DTOFs (f ECL ). coefficient, reduced scattering coefficient, and amplitude. The 2 cm data were analyzed through three rounds of fitting, with all three parameters allowed to vary freely in the first round. Thereafter, the mean amplitude for all wavelengths was computed, and the DTOFs were fit a second time with the amplitude fixed to this mean value, to obtain estimates of the wavelength-dependent absorption and reduced scattering coefficients. The wavelength-dependent scattering coefficients recovered from this round were then fit to a model of Mie Scattering (Eq. 1) to further reduce the noise in the estimated scattering coefficients. Note that in Eq. 1, a and b are the scattering amplitude and power, respectively, λ is the wavelength in nm, and µ s (λ) is the wavelength-dependent reduced scattering coefficient Thereafter, the DTOFs are fit a third time-with both the mean amplitude and reduced scattering values fixed to their values obtained from the previous round of fitting-to estimate the wavelengthdependent absorption coefficient (Figure 2).
The absorption coefficient depends on the tissue chromophores' concentrations (C i ) and their wavelength-dependent extinction coefficients (ε i ): This relationship can be used to recover the concentrations of oxyhemoglobin (HbO 2 ) and deoxyhemoglobin (Hb) from the wavelength-dependent absorption coefficients by solving the following system of linear equations (Eq. 3): where WF is the water fraction and is assumed to be 80% (Oros-Peusquens et al., 2019).

2.1.2.
Step 2: Rough estimation of the concentration of oxyhemoglobin and deoxyhemoglobin in the brain Step 2 establishes preliminary rough estimates of the concentrations of oxyhemoglobin and deoxyhemoglobin in the brain. One of the benefits of using trNIRS is that the tail of the DTOF is highly sensitive to deep tissue absorption (Diop and Lawrence, 2013). As such, the tail of the DTOFs were used to obtain a rough estimation of the brain absorption coefficient: Wavelength-dependent absorption (A) and reduced scattering (B) coefficient obtained following each round of fitting the 2 cm distance data assuming a homogenous medium.
where c is the speed of light in a vacuum, and n is the refractive index of the tissue (assumed to be equal to 1.4). The tail of the DTOFs were fit within the region between 5 and 1% of the max after the peak. All wavelengths from the long sourcedetector distance (3 cm) were individually analyzed using Eq. 4. The absorption coefficients were used to compute the rough initial estimates of oxyhemoglobin and deoxyhemoglobin concentrations in the brain, similar to Eq. 3. These concentrations will be later used to constrain the fitting of the 3 cm data with the two-layer analytical model (Step 4).

Step 3: Estimation of the reduced scattering coefficient
Step 3 analyzes the DTOFs from the 3 cm source-detector distance, using the process outlined in Step 1, to estimate the mean reduced scattering coefficient of the combined two-layer (i.e., both the ECL and the brain). Note that we are assuming homogenous tissue scattering, instead of layer-specific scattering coefficients, as this reduces both the duration and complexity of the fitting.

Step 4: Refining the extracerebral concentrations
In Step 4, the estimates of oxyhemoglobin and deoxyhemoglobin concentrations in the ECL are refined by simultaneously fitting the DTOFs of all the wavelengths with an analytical model of light transport in a two-layer diffusive medium (Kienle et al., 1998). The fitting parameters for this step are the concentrations of oxyhemoglobin and deoxyhemoglobin in both layers and the ECL thickness. The fitting algorithm was implemented in MATLAB, using the function fminsearchbnd to minimize the difference between the analytical model and the DTOFs of the full spectrum.
In this step, the concentrations of oxyhemoglobin and deoxyhemoglobin in the brain are fixed to the values obtained in Step 2, ECL concentrations estimated in Step 1 are used as the initial guess and allowed to vary ± 50%, and the ECL thickness is assigned an initial value equal to the average from the segmented MRI and allowed to vary ± 10%. The reduced scattering coefficients were fixed to the values obtained in Step 3. Finally, the amplitude is assigned by multiplying the one-layer amplitude from the homogenous fitting in Step 3 with the quotient of the one-layer solution divided by the two-layer solution, using the homogenous optical properties from Step 3. This corrects for any discrepancy in the amplitude of the two solutions. The fitting was conducted on the 3 cm source-detector distance DTOFs from 50% of the maximum on the leading edge to 5% after the peak. The refined ECL concentrations and thickness are then used in Step 5 to refine the brain concentrations.

Step 5: Refining the cerebral concentrations
In the final step of Phase 1, the ECL parameters recovered from Step 4 are fixed, and the concentrations of oxyhemoglobin and deoxyhemoglobin in the brain "layer" are allowed to vary between 0 and 80 µM, with an initial value of 30 µM. Note that the rough estimates of brain chromophores' concentrations obtained in Step 2 could be used as initial guess and allowed to vary ± 50%. However, because this step is very stable once the scattering coefficient, the ECL thickness, and the concentrations of the chromophores in the ECL are determined (Step 3 and 4), we relaxed the constraints for the concentrations of oxyhemoglobin and deoxyhemoglobin to illustrate the flexibility of the algorithm. The DTOFs obtained at the 3 cm source-detector distance for all wavelengths were simultaneously fit using the two-layer analytical model, from 5 to 1% of the maximum on the tail edge of the curves. This choice was guided by the fact late-arriving photons (i.e., the tail of the DTOF) are more sensitive to deep-lying tissue which, in this case, represents the brain layer. The estimated concentrations represent the baseline brain concentrations of oxyhemoglobin and deoxyhemoglobin, which are used in Phase 2 for rapid estimation of subsequent cerebral blood content and oxygenation.

Phase 2: Rapid estimation of cerebral blood content and oxygenation
Phase 2 uses the absorption coefficients estimated using the latephoton analysis (Eq. 4) on the 3 cm source-detector data to correct for the ECL contamination. This is accomplished by leveraging the absolute concentrations of the chromophores in both layers measured in Phase 1 (Eq. 5). Phase 2 has two steps: estimation of the ECL contribution to the absorption coefficient computed using latephoton analysis of the 3 cm SDD data (Step 6), and rapid estimation of chromophore concentrations in the brain (Step 7).
Step 7 can be repeated for all subsequent datasets (i.e., spectrum of DTOFs), enabling rapid neuromonitoring (∼0.18 s per analysis). Note that while Phase 1 is compatible with a spectral resolution of 10 nm, Phase 2 requires at least a spectral resolution of 2 nm between 680 and 930 nm. The increased spectral resolution is needed for the latephoton analysis due to the reduced signal-to-noise ratio (SNR) of late arriving photons.

Step 6: Estimation of the extracerebral contribution
Step 6 estimates the contribution of the ECL to the absorption coefficients recovered from the late-photon analysis of the DTOFs from the 3 cm SDD. This allows the fitting algorithm to account for the ECL contamination in subsequent hyperspectral trNIRS datasets, without the need to repeat the more computationally intensive twolayer analysis. More specifically, the absorption coefficients obtained from the late-photon analysis, µ a LP (t, λ), can be expressed as a weighted average of the absorption coefficients of the two layers: (Jacques, 2013).
where µ a Brain (T = 0, λ) and µ a ECL (T = 0, λ) are the baseline brain and ECL absorption coefficients computed by inputting the absolute concentrations estimated in Step 4 and 5 into Eq. 2, and assuming 80% water concentration. µ a LP (T = 0, λ) are the absorption coefficients recovered by applying the late-photon analysis to the baseline DTOFs from the 3 cm SDD. In a two-layer head model, the contribution of one layer can be expressed as a function of the other, such as f Brain = 1 − f ECL , where f Brain is the fraction of the absorption coefficient coming from the brain layer and f ECL is the fraction from the ECL. Thus, Eq. 5 can be rearranged to isolate the ECL fraction: f ECL (λ) = µ a LP (T = 0, λ)−µ a Brain (T= 0,λ) µ a ECL (T = 0,λ)−µ a Brain (T= 0,λ)

2.2.2.
Step 7: Rapid estimation of deep tissue chromophore concentration Step 7 allows rapid estimation of oxyhemoglobin and deoxyhemoglobin concentrations in the bottom layer (i.e., the brain). By combining Eqs. 5 and 6, the time-dependent absorption coefficient in the brain, µ a Brain (T, λ), can be expressed as a function of the baseline absorption coefficients in each layer, f ECL , and the baseline late-photon absorption coefficient: Furthermore, using Eq. 2 and assuming a water concentration of 80%, the wavelength-dependent absorption coefficient can be used to compute the concentrations of oxyhemoglobin and deoxyhemoglobin in the brain. Cerebral oxygen saturation (SO 2 ) and total cerebral hemoglobin (HbT) are then computed using Eqs. 8 and 9, respectively: HbT = C HbO 2 +C Hb (9)

Validation
Validation was conducted using in silico data generated with Monte-Carlo Extreme (MCX) in an adult head model at various skin and brain oxygen saturations (Fang and Boas, 2009;Yan and Fang, 2020). Using in silico data allows for direct comparison of the results of the algorithm with the "ground truth" inputted parameters while maintaining realistic geometry and optical properties. To better mimic experimental conditions, we segmented an MRI of an adult head into four tissue types: scalp, skull, cerebrospinal fluid (CSF), and brain using 3DSlicer (Kikinis et al., 2014). Brain and scalp oxygen saturations were varied independently from 40 to 80% and 50 to 70%, respectively, in 2% increments. The wide range of scalp oxygen saturations were simulated to investigate the confounding effects of changing scalp SO 2 on the accuracy of estimating brain SO 2 . 126 simulations were conducted for each brain-skin pair, corresponding to the wavelength range of 680 to 930 nm, in 2 nm increments. In total, we completed 8,316 simulations for this validation. The source was positioned on the right side of the head with the detectors placed 2 and 3 cm toward the forehead, as shown in Figure 3. Each simulation had a total of 3 billion photons with random seeds, ensuring realistic photon statistics with high SNR at both detectors. The optical properties of bone and CSF, the scattering coefficient, anisotropy factor, and refractive indices of the skin and brain, were obtained from literature (Firbank et al., 1993;Jacques, 2013). Absorption coefficients of skin and brain were computed using Eq. 2 for each oxygen saturation in the aforementioned range. Total hemoglobin in the skin and brain were set to 12.4 and 55 µMol, respectively, for all the simulations (Luttkus et al., 1995;Kanti et al., 2014;Auger et al., 2016;Farzam et al., 2017).

Results
Estimation of the baseline chromophore concentrations from a full spectrum of DTOFs (i.e., Phase 1) required 315 s of computation (CPU: Intel Core i7-6800K @ 3.4 GHz, using parallel computing Adult head model with source (red dot) and detectors (green dot). The source-detector distances were 2 and 3 cm. The tissues shown are skin (dark blue), skull (light blue), cerebrospinal fluid (orange), and brain (yellow). Due to GPU limitations, only the right upper octant of the head was simulated and is shown above-gray tissues not included in the simulations.

Phase 1: Baseline chromophore concentrations
The estimated absorption spectra and the expected wavelengthdependent absorption for four brain and three skin oxygen saturations are shown in Figure 4 to illustrate the qualitative similarity between the recovered and "ground truth" absorption spectra. Figure 4 shows that when cerebral oxygen saturation decreases, the error in the recovered absorption spectrum increases slightly; this is particularly obvious above 800 nm and is a result of errors in the recovered oxyhemoglobin concentration ( Figure 4A). Further, Figure 4B shows that the algorithm has excellent sensitivity to the brain as the ECL oxygen saturation has negligible impact on the accuracy of the estimated cerebral absorption coefficient.
The results of the quantitative analysis are shown in Figure 5A. The mean (± standard deviation) difference between the recovered brain oxygen saturations and the true values is 2.7 ± 2.5%, and the correlation coefficient between the true and recovered values is 0.99 (p < 0.0001). For the total hemoglobin, the mean (± standard deviation) difference between the recovered concentration and the simulated concentration of 55 µM is 2.8 ± 1.8%. The mean ECL thickness was estimated to be 12.4 mm, representing a mean (± standard deviation) difference of 3.0 ± 2.2% from the "ground truth." For comparison, Phase 1 was repeated with the ECL thickness fixed to its known value of 12 mm, which led to more accurate estimates of all the hemodynamic parameters ( Figure 5B). Notably, the mean (± standard deviation) difference between the recovered brain oxygen saturation was reduced to 1.5 ± 1.4%, with a Pearson's correlation coefficient of 0.99 (p < 0.0001), and the total hemoglobin was recovered with a mean (± standard deviation) difference of 1.7 ± 1.1%.
Correlation and Bland-Altman plots analysis were conducted for both when the ECL thickness was unknown and known (Figure 6). When ECL thickness was unknown (Figure 6A), the correlation plot had a sum of squared error of 2.3% and a Pearson R-value squared of 1, while the correlation plot for the known ECL thickness ( Figure 6C) had a sum of squared error of 1.3% and an R-value squared of 1. The Bland-Altman plot for the unknown ECL thickness case ( Figure 6B) had a coefficient of variation of 2.0%, while the plot for the known ECL thickness ( Figure 6D) had a coefficient of variation of 0.81%.

Phase 2: Rapid estimation of cerebral chromophore concentration
The estimated cerebral oxyhemoglobin, deoxyhemoglobin, and total hemoglobin from Phase 2 are shown in Figure 5C. The mean (± standard deviation) difference between the recovered and true brain oxygen saturation and total hemoglobin are 1.5 ± 1.5 and 3.1 ± 0.9%, respectively. Further, there is a strong agreement between the estimated and simulated values, with a Pearson's correlation coefficient of 0.99 (p < 0.0001). Figure 7 shows the correlation and Bland-Altman plots for the results of Phase 2. The correlation plot had a sum of squared error of 6.2% and a Pearson R-value squared of 1.00. The coefficient of variation from the Bland-Altman plot was 0.96%.

Discussion and conclusion
A hyperspectral trNIRS data analysis method that can quickly and reliably estimate cerebral concentrations of oxyhemoglobin and deoxyhemoglobin was developed and validated. The findings (A) Estimated absorption coefficient for the 40, 50, 60, and 70% brain oxygen saturations with skin oxygen saturation set to 70% versus the true (inputted) values. (B) Recovered versus true (inputted) absorption coefficient spectra for 50, 60, and 70% skin oxygen saturation at 70% brain oxygen saturation. Results from Phase 1 with unknown (A) and known (B) ECL thickness and Phase 2 (C): For every cerebral oxygen saturation, 11 scalp saturations (50 to 70 in 2% increments) were evaluated. The estimated oxyhemoglobin (red), deoxyhemoglobin (blue), and total hemoglobin (black) are plotted against the true cerebral SO2. The dotted lines are the "ground truth" while the solid lines are the estimated values. The error bars represent the standard deviation.
of the current report show that the algorithm can estimate adult cerebral oxygen saturation and total hemoglobin concentration with high accuracy. Although it takes 5.2 min to analyze a full trNIRS spectrum with Phase 1, this phase is only needed once per subject to establish the baseline concentrations of oxyhemoglobin and deoxyhemoglobin. Once concentrations are established, all subsequent trNIRS spectra can be analyzed with Phase 2, which takes only 0.18 s, to estimate cerebral hemoglobin concentration Correlation and Bland-Altman plots for the Phase 1 results with unknown (A,B, respectively; blue), and known (C,D, respectively; red) ECL thicknesses. The correlation plots show the mean estimated cerebral oxygen saturation for every true brain oxygen saturation (coloured squares), the line of best fit (solid black line), and the expected values (dotted black line). Correlation (A) and Bland-Altman (B) plots of Phase 2 results. The correlation plot shows the estimated mean cerebral oxygen saturation for every true brain oxygen saturation (blue squares), the line of best fit (solid black line), and the expected values (dotted black line). and oxygen saturation. The increased speed allows for real-time neuromonitoring while maintaining high accuracy.
A major challenge when using NIRS for adult neuromonitoring is the significant signal contamination by the ECL. Previous studies have reported that the ECL is responsible for 52-88% of the detected optical signal at 2-3 cm source-detector distances, which represents a serious challenge when estimating brain chromophore concentrations (Selb et al., 2014a;Milej et al., 2020a). To investigate Results of Phase 2 analysis when using a wide range of brain and scalp oxygen saturations as the baseline values for estimating f ECL . In the legend, (B) denotes the simulated cerebral oxygen saturation and S represents the simulated scalp oxygenation, both in percent. (A) The recovered oxyhemoglobin, (B) recovered deoxyhemoglobin, (C) recovered brain oxygen saturation, and (D) recovered total hemoglobin plotted versus the true cerebral oxygen saturation. The dotted black lines represent the "ground truth". the potential confounding effects of the ECL on the accuracy of the method, 11 scalp blood oxygenations were simulated for every brain oxygen saturation. As shown in Figure 5, the algorithm can estimate cerebral oxyhemoglobin and deoxyhemoglobin with high accuracy despite significant changes in scalp blood oxygen saturation. The reliability of the method under a variety of ECL conditions bodes well for its use in cardiac surgery and other high-risk procedures, since the brain is considered an index organ for oxygen supply to other organs; compromised cerebral oxygen supply may be indicative of larger systemic issues (Murkin, 2011).
To further explore the impact of different baseline concentrations on the accuracy of Phase 2, the method was tested using the extremes of our simulated scalp and brain oxygen saturations; cerebral oxygen saturations of 40 and 80%, and scalp oxygen saturations of 50 and 70%. The analysis revealed only minor deviations with the trends remaining unchanged, and all the tested scenarios resulted in the same Pearson's correlation coefficient of 0.99 for the recovered cerebral oxygen saturations (Figure 8). Further, a two-way ANOVA was conducted and did not reveal any statistically significant differences between the four datasets (p = 0.50 for the brain, and p = 0.49 for the scalp), showing that the fraction of optical contamination from the ECL does not significantly change with scalp or cerebral blood oxygen saturation for a given geometry within the tested range.
A major benefit of the method introduced in this report is that the exact ECL thickness does not need to be known precisely, as the fitting algorithm estimates (in Step 4) the ECL thickness within 10% of the initial guess. This could be beneficial during emergencies when there is not sufficient time for patients to receive an MRI or CT scan that would be used to estimate the ECL thickness. In such a scenario, the initial guess of the ECL thickness could be set to an average adult value as the algorithm will estimate it within 10% of the initial guess. Additionally, the relaxed requirement for the a priori knowledge of the ECL thickness reduces the impact of potential errors in its estimation from medical images on the accuracy of the method, as minor variations in ECL thickness are accounted for by the algorithm.
To improve robustness and reduce crosstalk, we assumed a cerebral water concentration of 80%. Since the water concentration in our simulations was 80% as well, we conducted additional analyses by assuming cerebral water concentrations of 70 and 90%. The analysis revealed that such error, which is ± 10% of the true water concentration, has a negligible impact on the accuracy of the recovered oxyhemoglobin and deoxyhemoglobin.
To reduce computational burden, Phase 1 uses a sparse spectrum (10 nm spectral resolution) instead of the dense spectrum (2 nm spectral resolution) used in Phase 2. The sparse sampling had negligible effects on accuracy; however, accuracy decreased when the spectral sampling was further reduced (i.e., more than 10 nm separation between consecutive wavelengths). Furthermore, the algorithm can be modified to recover other chromophores by increasing the number of chromophores in the fitting parameters and changing the wavelength region of interest to cover the spectral range of the target chromophores. Notably, future work will include adapting the algorithm to monitor the redox state of Cytochrome C Oxidase (CCO), a key biomarker of cerebral oxidative metabolism (Rajaram et al., 2018(Rajaram et al., , 2022Lange et al., 2019).
A potential limitation of this work is that it did not include the effects of the IRF in the Monte Carlo simulations. When analyzing trNIRS measurements, the IRF can be easily measured and its effects accounted for by convolving it with the analytical model. Note that for hyperspectral trNIRS measurements, the IRF of each wavelength must be measured as it has been shown that both the shape and temporal position of the IRF change with wavelength (Ioussoufovitch et al., 2021).
Furthermore, the initial rough estimates of the chromophores' concentrations in the ECL were obtained by assuming that 2 cm SDD mainly probed the ECL. To investigate the validity of this assumption we determined the partial pathlength of all the photons detected at the 2 cm SDD in the brain and ECL, for the four extreme cases of scalp and brain oxygenation (i.e., brain oxygenation at 40 and 80% and scalp oxygenation at 50 and 70%). The analysis revealed that more than 80% of the optical pathlength of the photons detected at the 2 cm SDD are in the ECL. Therefore, the optical properties obtained by analyzing the 2 cm SDD data with a semi-infinite homogeneous model are heavily weighted toward the ECL parameters and should provide a good initial guess of the ECL optical properties.
In contrast to Steps 4 and 5, wherein all the DTOFs were fit simultaneously, Steps 1, 2, and 3 were not spectrally constrained. The authors acknowledge that spectrally constrained approaches such as the method reported by D' Andrea et al. (2006) could be used in Steps 1, 2, and 3. However, given that these estimates were only used to obtain an "educated" initial guess of the concentration of light absorbers in the ECL, which were refined in Step 4, it is reasonable to assume that applying the spectrally constrained method would not significantly alter the accuracy of the algorithm.
Another potential limitation of this work is that the simulations were conducted on one octant of the full head due to GPU memory limitations. To assess the potential impact of this approach, we randomly picked 20 simulations and analyzed the photons' trajectories. We found that no detected photons reached the boundaries of the octant in the period of interest; thus, it is reasonable to treat the volume used in the simulations as an optically semiinfinite medium.
An additional limitation is that all validations were conducted in a single geometry. While the oxygen saturations in the scalp and brain were varied extensively, the medium was kept constant throughout. Nevertheless, we expect the algorithm to work with other ECL thicknesses since this parameter is estimated in Step 4 and thus does not need to be precisely known. Future work will include further validation in tissue-mimicking phantoms with various top layer thicknesses and in a pig model of an adult head before human applications.
In summary, this report introduces a multi-step hyperspectral trNIRS data analysis method that allows for accurate estimation of cerebral hemoglobin content and blood oxygenation in adults. Importantly, the approach does not require a priori knowledge of the ECL thickness, which increases robustness and usability. Further, the method can provide rapid estimates of cerebral oxygen saturation and total hemoglobin content once the baseline values are known, enabling real-time neuromonitoring. It is anticipated that this method would be valuable in a wide range of applications for continuous adult neuromonitoring, including cardiac surgery and intensive care.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
DC conducted the Monte Carlo simulations used for the validation and was the lead contributor to the development of the analysis method described. NL helped develop Phase 2 of the analysis method. SI provided feedback during the development of Phase 1. MD was the primary investigator for this project, providing insight and guidance throughout, as well as securing funding. All authors contributed to the article and approved the submitted version.